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We study the emergence of Boltzmann's law for the "single 
particle energy distribution" in a closed system of interacting 
classical spins. It is shown that for a large number of parti- 
cles Boltzmann's law may occur, even if the interaction is very 
strong. Specific attention is paid to classical analogs of the 
average shape of quantum eigenstates and "local density of 
states" , which are very important in quantum chaology. An- 
alytical predictions are then compared with numerical data. 

PACS numbers: 05.45.-a 



I. INTRODUCTION 

Chaotic properties of few degrees of freedom systems 
have attracted much attention during the last years. The 
knowledge and classification of chaotic dynamical sys- 
tems is now extremely accurate, if compared with the 
status of art at the beginning of the century. On the 
other hand, conventional statistical mechanics was born 
long before the chaology of dynamical systems. It is well 
known that, in contrast to chaotic dynamical systems 
with few degrees of freedom, for the onset of statistical 
equilibrium in systems with an infinitely large number 
of particles neither the non-linearity nor the interaction 
between particles is needed. Indeed, in such systems 
the thermodynamic limit (infinite number of particles, 
N — > oo ) leads to the statistical behavior of a system 
since any weak interaction (with an environment or be- 
tween particles) gives rise to mixing properties and onset 
of chaos. 

Although statistical mechanics has obtained undoubt- 
edly results in many different fields, fundamental ques- 
tions about statistical description of systems with fi- 
nite number of particles are still open Q. It is clear 
that in such systems the interparticlc interaction is cru- 
cial, however, its role in producing chaos or ergodicity is 
still not understood. Due to computing difficulties and 
lack of theoretical studies, many-body chaotic systems 
have been scarcely investigated and little is known about 
them. 

Recently, a novel approach to quantum isolated sys- 
tems of finite number of interacting particles has been 
developed |^-^]. The goal of this approach is a direct 
relation between the average shape of exact eigenstates 



(F— function), and the distribution of occupation num- 
bers n s of single-particle levels. This relation shows 
that there is no need to know exactly the eigenstates, 
instead, if these eigenstates are chaotic (random super- 
position of a very large number of components of basis 
states), the F— function absorbs the statistical effects of 
interaction between particles and determines the form of 
the n s — distribution. The approach has been mainly de- 
veloped for a model with completely random two-body 
interaction between a finite number of Fermi-particles, 
however, it can be also applied to dynamical models with 
chaotic behavior. 

One of few dynamical models studied with the use of 
this approach is a system of two interacting spins 
The most interesting result obtained numerically in Ref. 
|Q , is that the distribution of occupation numbers (even 
for two interacting particles) in the chaotic region can be 
described in the form of the standard Bose-Einstcin (BE) 
distribution (only symmetric states have been considered 
which results in the BE statistics). It was shown that the 
statistical effects of the interaction lead to an increase 
of temperature which is related to the BE distribution. 
Also, it was discovered that the canonical distribution 
is recovered if one randomize the non-zero elements of 
the interaction V keeping the dynamical constrains of 
the model. This means that random interaction plays 
the role of a heat bath and allows to use statistical and 
thermodynamical description for isolated systems. 

In order to extend the approach of Refs. (|-||] to clas- 
sical dynamical systems with a large number of particles, 
we introduce below the model of A~ interacting spins in 
one-dimension. Due to the well-defined classical limit, 
it is of the first interest to explore similar problems in 
the classical counterpart. The problem of a quantum- 
classical correspondence for chaotic systems with a large 
number of interacting particles seems to us extremely 
important in view of many physical applications. In this 
paper we mainly concentrate on the occurrence of Boltz- 
mann's law and on the classical counterparts of quantum 
local density of states and shape of eigenfunctions. We 
deserve for future work [pj the investigation of the anal- 
ogous Fermi-Dirac or Bose-Einstein statistics in a quan- 
tized version of such a model. 

Our investigation is complementary to the approach 
based on the so-called dynamical temperature 1 10-1^], in- 
troduced for the study of the statistical properties of few 
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degrees of freedom classical models with chaotic behav- 
ior. Specifically, our interest is in the notion of statistical 
temperature associated with Boltzmann's distribution (if 
any) for single particle energy. The relation between sta- 
tistical and dynamical temperatures represents an open 
and interesting question. 



II. THE MODEL OF INTERACTING SPINS 

Our purpose is to investigate a few-body dynamical 
system from the statistical mechanics point of view. Fol- 
lowing previous works on two-particle spin problem ^[t| , 
we consider the Hamiltonian: 
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which is known as the 1-d Heisenberg model in a mag- 
netic field, see for instance |U5). Here Si are spin vectors, 
and, for sake of definiteness we take periodic boundary 
conditions, S\ = <S/v+i- The classical version is a solvable 
model, from the statistical point of view (namely when 
the number of particles N — > oo, see ]l4[]). 

Instead, we are here interested in a dynamical ap- 
proach, in particular, for a small number of particles, 
when the usual statistical approach is at least question- 
able. In what follows, we shall consider a more simple 
version of the classical Heisenberg model, which is de- 
scribed by the following Hamiltonian : 



H = H + V = BYZi S ! + J ^li S ? S i+i ( 2 ) 

where N > 2 is the number of spins in the chain. This 
we do, in order to simplify analytical calculations. The 
equations of motion can be written in the usual way: 



(3) 



where {, } are the Poisson brackets (see |15|,|16[ for similar 
dynamical models). Constants of motion are the energy 
E and the magnitude of the angular momenta \Si\ = s 
(the latter assumed to be the same for each spin) . With- 
out loss of generality, we can put s = 1. 

Contrary to the common view point (spin Hamiltonian 
plus a magnetic field as a perturbation) we consider the 
"magnetic" part as an unperturbed Hamiltonian. Indeed, 
we are interested in models which can be generally ex- 
pressed as a sum of single particle Hamiltonians, this fea- 
ture not being shared by the perturbation V. Therefore, 
the Heisenberg kernel (J SiSi+i) wm be considered as 
a perturbation between nearest-neighbor spins. 

The unperturbed Hamiltonian is integrable and the so- 
lution of the unperturbed equations of motion can be 
written down at glance. It should be pointed out that 
the perturbation itself is not chaotic. Indeed, it can be 



verified numerically that two close trajectories diverge 
only linearly in time when B = 0. This means that the 
maximal Lyapunov exponent is zero for any choice of the 
initial conditions. On the other side the same Lyapunov 
analysis shows that the total Hamiltonian, for generic val- 
ues of the coupling J and for energy values in the middle 
of the band [E ~ 0) is chaotic. Maximal Lyapunov expo- 
nents, as a function of time for different J and B values 
and fixed energy E , are shown in Fig.|l|. 
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FIG. 1. Maximal Lyapunov exponents as a function of 
time for E — — 2, J = 1 and different B = 0.5, 1 as indicated 
in the picture. The system has N = 4 spins. 

As one can see, the Lyapunov exponents approach zero 
when B decreases, or J increases. On the other side, for 
fixed J and B different from zero, the region where the 
maximal Lyapunov exponent approaches zero is close to 
the edge of the energy spectrum (see Fig.||). Maximal 
Lyapunov exponents have been calculated using the stan- 
dard recipe Jl7| . 
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FIG. 2. a) Maximal Lyapunov exponents as a function of 
time for J = 1, B = 1 and different energy E as indicated in 
the text. Simulations have been done with 4 spins. 



Besides integrability, 
Hq, in the limit N — > < 



the unperturbed Hamiltonian 
o , has good statistical proper- 
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ties which we are going to exploit in the next Section. It 
is, indeed, interesting to compare our results with those 
obtained with the standard statistical approach. 



III. IDEAL GAS OF SPINS 

This particular choice of the model, being extremely 
simple, allows an analytical treatment of the case N — ► oo 
in both unperturbed and perturbed case. Let us analyze 
first the unperturbed Hamiltonian Hq. This should be 
thought as a model for a very weakly interacting sys- 
tem. In what follows, without loss of generality, we as- 
sume B = 1. The microcanonical ensemble represents 
the most natural way to analyze an isolated system of 
N spins. One question is to find the thermodynamical 
temperature in dependence on the energy E . Another 
is to study the single particle energy distribution and 
the conditions under which it can be assimilated to the 
Boltzmann distribution. 

The temperature T — /3 _1 (ks = 1 hereafter) can be 
defined via the microcanonical ensemble: 



dS 



(4) 



where S is the entropy. For a sufficiently large number 
of particles this can be in turn defined as S = log p(E). 
Here p(E) can be defined through the phase space vol- 
ume If the motion of TV spins is ergodic, each of 
them covers uniformly the unit 3-d sphere. Therefore, 
each component of 5* has a uniform probability density 
function in the interval [—1,11 1191, that is : 



p(h) 



dPjS[ < h) 
dh 



1/2 when - 1 < h < 1 
elsewhere 



(5) 



and the same for the other components Sf and Sf . Here 
the quantity P(x < a) gives the probability that the 
continuous random variable x gets values less than a. In 
the same way, the density of states can be evaluated as 
a probability : 



dP(Eli S[ < Ep) 
dE 



Po(E ,N) 



(0) 



The distribution of the sum of N independent random 
variables can be obtained using the central limit theorem 
(when N — > oo), 



Po(E ,N) 
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As a result, from Eq. (|4j) one gets an unusual microcanon- 
ical relation for the energy vs. inverse temperature fll4|: 
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In principal, the energy Eq ranges from _E m i n < to 
B max = —E mnl > 0, therefore, in this model negative 
temperatures are possible too. One should note that 
typical physical systems have density of states increas- 
ing with energy, thus giving positive temperature. For 
this reason, in what follows, we consider only positive 
temperatures which correspond to the left (part of the 
symmetric) energy spectrum, Eq < 0. 

It is important to notice that the density of states can 
also be obtained for finite system of particles without 
invoking central limit theorems. It can be done by a 
direct integration over the phase space volume, 



p (E ,N) = C dSl p(Sf) . . . dS% p(S* N )x 
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°[J\ dx \ e lXx ] 
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(10) 

where M is the integer part of (Eq + N)/2. 

Let us now concentrate on the "single-particle energy 
distribution" (SPE- distribution) defined as : 



dP(S!<h\j:? =1 S?=E ) 
dh 



n (h\E ) = N 



(11) 



where P(x < a\y = b) is the conditional probability. 

This quantity is the classical analog of the quantum 
distribution of occupation numbers, which gives the num- 
ber of particles occupying a single-particle level with en- 
ergy Eq. Correspondingly, the above classical distribu- 
tion determines the probability that any of iV spins has 
the energy Eq. From this definition, it is clear that the 
SPE-distribution ( |TT1) is normalized to the total number 
of particles, and it defines the mean energy of all parti- 
cles, or, the same, the mean total energy of the system. 
In this way, one can treat an isolated model in the same 
way as the model in contact with a heat bath (see details 
inRefs. |, |). 

It is well known that for N — > oo the single particle 
energy should be distributed according to Boltzmann's 
law. Under suitable conditions this holds true in this 
model, too. Indeed, one gets 
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n (/i|£ ) = TV pot^Af) 



N 



TV 



Pa(E ,N) 



(12) 



_ at po(/t,l)po(g -/i,jV-l) 
- JV Po(E ,N) 

where the last equality is due to the independence of S"f . 
The quantity p(Sf — h, X^i" 1 &f — Eo) defines the joint 
probability density function. 

Substituting the expression (Q) into the above relation, 
it is easy to obtain (in the limit of a very large N ^> 1) : 



n Q (h\E ) ~ exp 
where E < 0. 
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FIG. 3. Single particle energy distribution for TV = 10, 
and J — 0. Full lines represent the exact distribution ( |l2) o b- 
tained from the exact unperturbed density of states (Eq.(lCj)). 
Dashed lines are the correspondent infinite approximations 
(^) with a proper normalization. Upper curves are for 
Eq = —6, lower ones for Eo = —2. 

Obviously, a correspondent expression for any num- 
ber of particles can be obtained, using the exact value of 
the unperturbed density of states Eq. (|lfj|) . In Fig.|| we 
show the infinite (Gaussian) approximation (TV — » oo) 
and the exact one, computed for two different energy 
values. As one can see while in the middle of the spec- 
trum the distributions are almost the same, close to the 
band edges they are remarkably different. The Gaus- 
sian, or infinite approximation, works very well even for 
a small number of particles as soon as the energy Eo is 
not close to the edges of the band (the interval [—TV, TV] 
for non- interacting spins). Indeed, the Gaussian density 
of states has infinitely long tails, while the exact one is 
sharply defined within [—TV, TV]. On the other hand while 
the support of the density function scales as TV, its vari- 
ance depends on \/TV. It follows that the region close to 



the edges becomes less and less important as TV becomes 
large. 

Let us briefly comment Eq.(fi"3|). Although it has been 
obtained for a large number of particles, nevertheless it 
shows that, generally speaking, the single particle en- 
ergy distribution for finite weak-interacting systems is 
not necessarily described by an exponential law. One can 
see that the latter occurs in a very strong limit TV — > oo, 
provided that \E \ 3> \h\ , 



n (h\E ) ~ exp(-/3 /i) 
where the temperature To = P^ 1 is defined via 
3£o 
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TV - 1 



P mc + 0(1/TV) 



(14) 



(15) 



and Eq < 0. One should stress that, apart from the limit 
TV — > oo , the exponential distribution arises when the 
second condition Eo 3> \h\ ~ 1 is fulfilled. Physically 
this means that the total energy must be larger than 
the typical single particle energy, a condition naturally 
satisfied for a thermodynamical system. 

Another important relation can be obtained from 
Eq. ( |r2[). S pecifically, taking the derivatives of both sides 
of Eq. (|14|) over h , and exchanging the derivative in the 
r.h.s., one gets : 



~~ -jr log n (h\E ) 
an 



dE 



logp (E -h,N-l) (16) 



One can see that in the limit TV — > oo and \Eo\ 3> 
the r.h.s coincides with the microcanonical definition of 
temperature while the l.h.s. shows that the only possible 
exponential distribution should have the microcanonical 
temperature. 

We have thus found that in the non-interacting sys- 
tem, when the number of particles becomes sufficiently 
large and the total energy of the system is larger than the 
typical single particle energy, the microcanonical temper- 
ature @ defined by the density of states, coincides with 
the statistical temperature defined directly from the SPE- 
distribution ( |IT| ) which is the standard Boltzmann law. 
One can then assume that if the interaction is sufficiently 
small, in order not to change dramatically the previous 
results, but enough to produce ergodicity (from the equa- 
tions of motion) , Boltzmann's law again results from tak- 
ing a sufficiently large but finite number of particles. 

Let us stress that in order to get these results, the mo- 
tion has been assumed to be fully ergodic on the unit 
3-d sphere. Rigorously speaking, this can be true, from 
the dynamical point of view, only in the presence of in- 
teraction J ^ 0. Indeed for J = the motion is foliated 
onto regular tori and no ergodicity at all is allowed. So, 
this ideal spin model has been worked out following the 
traditional statistical mechanics picture, where the in- 
teraction is assumed to be sufficiently weak in order to 
apply non interacting results, but sufficiently strong in 
order to get ergodicity. 
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IV. MANY INTERACTING SPINS 
A. Density of states 

1. Infinitely strong interaction 

Before switching to the interacting case, let us compute 
the density of states in the presence of an infinitely strong 
interaction (namely putting formally B = 0). In this case 
the density of states can be written as follows, 



dP(V<E) _P(JEl 1 S?Sf +1 <E) 



Pv(E,N) = 



dE 



dE 



(17) 



The central limit theorem can also be applied in this case, 
keeping in mind that the normalized probability distribu- 
tion of the product z = xy of two uniformly distributed 
random variables x, y in the interval (—1, 1) is given by 



p v (E/J,l)=-(l/2J) log(|£/J|) 



(18) 



so that the variance is (z 2 ) — 1/9. Eq.(|l7|) then becomes 
(in the large TV-limit) : 



Pv(E,N) 



exp 



E 2 
2a 2 v 



with 



NJ 2 
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(19) 



(20) 



Note, that even in this case an explicit integral equation 
can be obtained for finite systems. Following the same 
steps as in Eq.([To|), one gets : 



Pv(E,N) = 



-\- CO 



-i\E 



s(X) 



A 



iV 



(21) 



where s(A) is the sine- integral function defined by 



s(X) 



dx 



sin x 



2. General case 

In the same way the density of states in the presence 
of both the interaction and the "kinetic" term can be 
computed assuming Hq and V are independent. Indeed, 
let us define : 



P(H + V < E) = 



dE' p(E') 



(22) 



On the other hand, if Hq and V are independent, their 
joint probability density function can be written as 



Ph ,v(Eq,E) = p (E )p v (E) 



(23) 



The probability density function of the sum of two inde- 
pendent random variables is thus given by [po[ : 



/+oo 
dE'p (E-E')p v (E') 
-CO 



(24) 



Substituting Eqs.(|7|,|17|), and performing simple Gaussian 
integrals one gets : 



p(E,N) 



where the variance 



1 
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exp 
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a 2 = 



N 
3^ 



1 + J 2 /3 



(25) 



(26) 



should be compared with Eq.(|20j). 

In Eq.(|2^) there are two different approximations. The 
first one is the Gaussian form for the density of states, 
which is valid when N — > oo. The second is the indepen- 
dence of the (random) terms H and V which, of course, 
cannot be true in general. For instance, a configuration 
with all spins aligned along the z-axis, for which Eq = N 
implies V = 0. However, if the motion is ergodic and 
the energy is not too close to the band edges, the second 
assumption can be considered as a good approximation. 



Classical analogs of quantum eigenstates and the 
LDOS 



Another important relation can be obtained linking the 
two quantities, shape of 11 eigenf unctions" (SE) and local 
density of states (LDOS), whose concepts were motivates 
by quantum mechanics (see for example, P,^l|,p|j6]| ) . The 
latter quantity, LDOS, also known in nuclear physics as 
strength function, is very important when describing the 
spread of the energy initially concentrated in a specific 
unperturbed state. The classical analogs of these func- 
tions have been introduced in Ref. |22| and recently ap- 
plied to dynamical models in Refs. [p|]2l|. The very point 
is that in the limit of N — > oo these two quantities can 
be explicitly found and a relation between them can be 
established. 

In the case of ergodic motion the classical analog of 
the shape of eigenstates (SE) can be defined as 



W E (E ) 



dP{H <E \H a + V = E) 
dEn 



(27) 



Correspondingly, the classical analog of the LDOS is 



WE (E) 



dP(H + V< E\H = Eg) 

dE 



(28) 
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It is very important that both the SE and LDOS can be 
computed more efficiently p|l using the equations of mo- 
tion. For the SE one has to choose a chaotic trajectory at 
some fixed energy E, compute the Ho(t) trajectory and 
sample the values of the unperturbed Hamiltonian Hq 
along this trajectory at some fixed time intervals. This 
procedure gives us the ergodic distribution inside the en- 
ergy shell constructed by a projection of the phase space 
of H onto H (see details in Pj2l|). In the same way, 
the classical LDOS can be numerically computed tak- 
ing a bunch of (regular) trajectories of the unperturbed 
Hamiltonian Hq = Eq and computing the correspondent 
spread of H(t) along these unperturbed trajectories. The 
sample of the values of the total Hamiltonian H (t) taken 
at given intervals of time results in the classical LDOS. 
Let us stress that even in the case of ergodic motion of the 
total Hamiltonian H , when only one single trajectory is 
needed in order to get the SE, an ensemble of trajectories 
of Hq is necessary in order to get a reliable result for the 
LDOS. This is a consequence of the integrability of Hq. 

Alternatively one can choose, as indicated by the defi- 
nitions (|27j) , (p8|) , many different initial conditions on the 
same energy surface and sum over them. It is clear that 
the two procedures (dynamical and taking the average 
over the phase volume) should give the same result in 
the case of ergodic motion. 

An important identity can be easily proven from the 
relations (|27|) , (p8|) . Let us write, 

P{Ho + V<E,Ho<E ) = 

= j E l dE' Q P(H + V < E, Ho = = (29) 

dE * p ( H ° + E \ H ° = E 'M E 'o) 

from which one can get the classical SE, 

we (E)pq(E q 
In the same way one can obtain the classical LDOS, 

W E (E ) P (E) = d2p(H °f E l d T V ^ E) (31) 

From above, the following relation emerges, 

w Eq (E) Po (E ) = W E (E )p(E) (32) 

Let us stress that the previous identity does not de- 
pends on the Gaussian approximation and it takes into 
account dynamical correlations too. This simple relation 
between classical SE and LDOS is very important in dif- 
ferent applications. Remarkably, only unperturbed and 
perturbed densities of states, po(E) and p(E), get into 
this relation. This allows to relate the shape of eigen- 
states to that of the LDOS in the corresponding quan- 
tum model, in a deep semiclassical region. In fact, the 



d 2 P(H a +V<E,H Q <E Q ) / 3 qx 



knowledge of these two functions leads to a novel semi- 
classical approach according to which it is easy to detect 
quantum effects of localization, see details in . 

It should be stressed that classical SE and LDOS, in 
essence, are ergodic measures for energy shells defined 
by the projection of Hq onto H (and vice versa) in the 
energy representation. 

Functions W E (Eo) and w Eo (E) can be considered as 
kernel operators transforming unperturbed quantities in 
total ones and vice versa. For instance, by integrating 
Eq.(p|) one has : 



p(E) 



and the converse, 



Po(Eq) 



dE w Eo (E)po(E ) 



dE W E {E ) p(E) 



(33) 



(34) 



In a certain way w and W can be considered one the 
inverse of the other. It is also easy to check that when 
J -> then We(E ) = w Eo {E)) = 8{E - E ). 

Assuming the Gaussian approximation for the densi- 
ties of states, and using Eq.(|32|) we can easily get an 
analytical expression for the classical SE : 



We(E ) = 



1 



where 



2 



exp 



NJ 2 



{Eq — E c 
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(35) 



(36) 



is the variance, and E c = vE is the center of the SE. 

In order to obtain the LDOS distribution, one can use 
the relation (Ell) from which one gets, 



we (E) 



1 



where 



exp 



NJ 



(Eq - Ef 
2al 



(37) 



(38) 



It is interesting to note that the above distribution, in 
the Gaussian approximation, coincides with Eq.(|l9[), 



weo(E) =pv{E-E ). 



(39) 



This relation can be also obtained independently, using 
the assumption of the independence of Ho from V. 

One can see that although the two Gaussians are dif- 
ferent, for a weak interaction J « 1, and then v ~ 1, 
they appear to be close one to each other. This fact is 
of a general nature and occurs in other models, see Refs. 

On the other hand, this result shows that, strictly 
speaking, even in the case of ergodic motion, one should 
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not expect SE and LDOS to be the same. As was found 
above, the relation between the SE and LDOS is given 
by Eq. (|32|) . One can obtain a very useful relation for 
the variances of the SE and LDOS (valid in the Gaussian 
approximation only), 



(40) 



We would like to stress again that both SE and LDOS 
have a proper meaning only in the case of ergodic motion. 
In other cases they strongly depend on initial conditions 
and on the energy, thus not allowing to use any statistical 
approach. 



C. Distribution of single particle energies and 
different temperatures 

For TV 3> 1, the microcanonical temperature T mc = 



(3 1 can be defined from the total density of states (25) 



_ aiog(p) _ 3E 



(41) 



with v = 1/(1 + J 2 /3). In fact, Eq.(|l]) determines the 
thermodynamical temperature since it establishes the re- 
lation between the temperature and total energy of a sys- 
tem. Now, we are going to find the statistical temperature 
associated with the SPE-distribution. First of all, let us 
notice that, if the variables Ho and V are independent 
the following approximate relation can be obtained : 



n(h\E) = ^<H\H + V=E) = 
= J dE Q n {h\E ) W E (E Q ) 



(42) 

Substituting Eqs.(|T3|,|35|), in Eq.(|2j) one gets, 

n(h\E) ~ exp (-53^ - + (43) 



where E c = vE. Note that when J = 0, then v — 1 and 
Eq.( |l3|) is recovered. 

One can see that in the limit \E\ ^S> \h\ ~ 1 and N ^> v 
the SPE-distribution can be approximated by the expo- 
nential dependence: 



n(h\E) 



where 



(3{E) 



3E 

' AT 



-l3(E)h 



V = (3„ 



(44) 



(45) 



The last equality, claiming that the statistical temper- 
ature coincides with the thermodynamical one, is by no 
means trivial. From one side, the presence of strong in- 
teraction suggests some statistical equilibrium property 
thus leading to the microcanonical predictions. On the 



other side we do not know if the single particle descrip- 
tion is valid since the interaction is strong. Even in the 
case when mean field approach would be possible (in such 
a way that part of the interaction becomes actually a sin- 
gle particle energy), the single-particle distribution would 
depend on the mean field and there are no reasons to ex- 
pect, a priori, that the Boltzmann distribution occurs 
with the same temperature as given by the microcanon- 
ical ensemble. 

Another way to find the statistical temperature which 
corresponds to the SPE-distribution, is as follows. First, 
we note that in the presence of interaction one has 



dh h n{h\E) ^E^E 



(46) 



that is different from the non-interacting case, for the 
latter the obvious relation holds, 



J dhh n {h\E ) = E a 



(47) 



Note that the SPE-distribution in Eq.(f46|) depends on 
the interaction, however, the "mean value" of the energy 
of all particles does not correspond to the energy of a 
system. This very fact allows to relate statistical effects 
of the interaction to an increase of the total energy (in 
application to quantum systems this approach was con- 
sidered in f§§|f). 

Specifically, in order to find the SPE-distribution 
n(h\E) from Eq. (J46|) , one needs to know the renormalized 
energy E. One can see that the Boltzmann distribution 
gives a correct result for an isolated system, if we take 
into account the shift of the energy = E — E which 
is due to the interaction between particles. In a sense, 
the (statistical) effects of the interaction are absorbed by 
the increase of the energy of a system, compared to the 
case of non-interacting particles. Therefore, pseudoran- 
dom interaction may be treated as an internal heat bath, 
thus giving rise to a statistical equilibrium. 

The shift A e can be found by assuming the absence of 
correlations between Hq and V in a way described above. 
One can show that 



E 



dE E W E (E Q ) = E c 



(48) 



where E c — vE is the center of the classical SE, and 
the last equality occurs in a strong limit N ^> 1 using 
Eq. (|35"|) . Therefore, the shift is given by the following 
relation, 



cr 2 J 2 
E-E = — ^E=-^-E 



(49) 



It is important that the shift is defined by the width 
a w of the LDOS and by the width of the unperturbed 
density. One should note that in quantum representation, 
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the variance is denned by the sum of squared off- 
diagonal matrix elements. Therefore, this energy shift 
can be found without diagonalization of huge quantum 
Hamiltonian matrices. Using this shift, one can find the 
temperature from Eq.(|4^) by assuming the Boltzmann 
dependence for n(h). 

One more way to obtain the temperature is via the un- 
perturbed density of states evaluated at the renormalized 
energy E 



dlogpo 



dE 



E„ = E 



3E 



Using the relation between E and E, one gets, 

R -JA-^ - R 
Pr - N - N V - p mc 



(50) 



(51) 



Once again, we should remind that apart from the limit 
N ^S> 1, the above relation is valid when neglecting the 
correlations between Hq and V . However, the smaller is 
the number of particles, the larger is the energy region 
where (dynamical) correlations will be important. 

It is now interesting to find the increase of tempera- 
ture AT = T — T due to the interaction, in comparison 
with the temperature T — fi^ 1 of the system with non- 
interacting spins, 



AT 

tT 



3 



(52) 



One can see that the relative increase of temperature is 
given in terms of the variance of the LDOS and the width 
of the unperturbed density only (see also Eq.([f9|). This 
fact seems to be generic, see discussion in 

More accurate relation for the definition of tempera- 
ture via the SPE-distribution, without the assumption 
of the Gaussian form of the density of states, can be ob- 
tained as follows. In close analogy to Eq. (|l6j) one can 
write, 



n(e\E) 



dP{Sl<t\H=E) 



1 dP(SZ<e,H=E) 
p(E,N) de 



(53) 



On the other hand, defining Xi 
it follows : 



Sf , and yi — J Sf Sf^_-^ 



P(xi <e,H = E) = 

= floe de'P(x 1 = e',x 1 +y 1 + J2i=2 Xi + Vi=E) = 

= de'P{ Xl = e>, Vl + Ejl 2 ^ + Vi = E - e') = 

= J dhf de'P{x 1 = e', y 1 = h, Ya= 2 Xi + Vi = E - e' - h) 

(54) 



Owing the independence of the 3 random variables, the 
joint probability density can be factorized and one gets : 



(55) 



n(e\E)p(E,N) = 

= po(e,l)/ dh p v (h,l)p(E-e-h,N-l) 
Defining now the "interacting" density of states 

Z e (E) = J dh Pv (h, l)p(E-e-h,N-l) (56) 
one gets a relation similar to Eq. (|l^ ) : 



-^logn(e| J B)= ^\og Z e (E) 



(57) 



Eq. ( p7\ ) states that if correlations can be neglected then 
a renormalized density of states must be introduced in or- 
der to get the "correct" temperature as obtained from the 
single particle energy distribution. On the other hand, it 
is clear that in the limit N — + oo, and \E\ ^> 1, we have 
Z e (E) ~ p(E) and the usual definition is recovered. It is 
now interesting to apply our estimates to a system with 
not very large number of spins. 



4.5 
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dynamical 




Gaussian 




x dynamical 




Gaussian 
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FIG. 4. Single particle energy distributions for N = 100 
and large interaction J = 1, for different energy values. One 
single trajectory has been iterated for a given energy, up to a 
time t — 10 6 . Numerical data should be compared with the 
analytical expression Ji3| ) obtained from the Gaussian approx- 
imation for the density of states. Circles stand for E = —19, 
crosses for E — — 7. 



V. NUMERICAL DATA 

A. Large number of spins 

Let us first consider the model with N — 100 spins. 
On one hand, this situation is far from the thermody- 
namic limit, on the other hand, the number of particles 
is quite large and one can expect a good correspondence 
with analytical results obtained in previous Section. In 
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Fig.|| two SPE-distributions are plotted for different val- 
ues of the total energy E. In both cases, the trajecto- 
ries have been found to be chaotic with positive maximal 
Lyapunov exponents. For comparison, the analytical ex- 
pression (see Eq.|43|), obtained in the Gaussian approxi- 
mation, is shown. As one can see, the agreement is fairly 
good. 




-25 -20 -15 -10 -5 



E 

FIG. 5. Inverse temperatures versus total energy for 
N — 100 and J — 1. Crosses are the extracted j3fu from 
the best fit of the SPE-distributions n(h\E) to the exponen- 
tial dependence. Typical examples of these distributions are 
given in Fig^. Dashed line is the microcanonical relation 
( pl| ) found numerically from the density of states p{E) . Full 
line is the theoretical prediction, obtained from Eq.d57|) with 
e = — 1. Dotted line is the Gaussian approximation (see in 
the text). 



of \nn(h\E) at the bottom of the energy spectrum, 

_d\nn{h\E) _ 3E C 3vh _ 3v(E - h) 
fit ~ dh ~ N-v ~ N - v ~ N-v 

(58) 

with h —* — 1. The corresponding result is shown in 
Fig.|| by the dotted line. One can see that this depen- 
dence practically coincides both with the fit to actual 
distribution n(h\E) (crosses), and with the temperature 
determined by Eq.(|56|) (full line). This means that these 
two approximations correctly take into account both the 
dynamical correlations and the finite number of particles. 

We can conclude that, in spite of the relatively small 
number of particles, if compared with the usual ther- 
modynamical systems of 10 23 particles, our model of 
N = fOO interacting spins can be approximately de- 
scribed by standard statistical approach which ignores 
dynamical correlations between particles. Indeed, the 
difference between microcanonical temperature and the 
approximate temperature found by the fit of the SPE- 
distribution close to the edge of the energy spectrum, is 
quite small and may neglected in some cases. However, 
even for a relatively large number of spins, N — 100 ^> 1, 
a clear influence of dynamical correlations and finite num- 
ber of particles remains. 

One should, however, stress that in spite of a clear 
manifestation of the influence of dynamical correlations 
and finite number of particles, the SPE-distribution 
n(h\E) for low energies can be effectively described by 
the standard Boltzmann distribution, however, with a 
renormalized temperature. 



On the other hand, a direct comparison of the micro- 
canonical temperature (3 mc with the approximated tem- 
perature, obtained by fitting the SPE-distribution by an 
exponent, gives different answers, see FigJ^. 

In this figure crosses represent the fitted inverse tem- 
peratures as a function of the total energy, while the 
dashed curve is the microcanonical relation between the 
inverse temperature and the energy. In order to smooth 
fluctuations, the derivative has been calculated perform- 
ing local averages in small energy windows. A small, but 
systematic, difference between the statistical and micro- 
canonical temperatures is clearly seen. 

A much better agreement can be achieved by using 
directly, (instead of the total density of states) the renor- 
malized density Eq. ( |56| ) . The corresponding result is pre- 
sented by the full line. One can see quite good correspon- 
dence to numerical data for the statistical temperature. 
This result means that Eq. (p&j is more accurate than the 
usual statistical definition (El). 

One can also introduce the effective (statistical) tem- 
perature by making use of Eq. ( |l3| ) , keeping all terms in 
the exponent. The temperature can be defined as a slope 



B. Chaos versus ergodicity 

It is reasonable to think that a statistically stable dis- 
tribution would require a certain degree of chaoticity. 
But chaos itself, as indicated, for instance, by the posi- 
tivity of the maximal Lyapunov exponent, is not enough 
in order to get equipartition among different degrees of 
freedom (see general discussion of this very important 
problem in |23]). We have found such cases for a small 
number of particles and energy close to the center of the 
spectrum. 

It is instructive to work out a specific example. Let us 
consider a system of N = 4 spins. For J = B = 1 and en- 
ergy E* = —1.9745 the single particle energy distribution 
is very different from the phase average distribution (ob- 
tained with many different initial conditions in a small 
energy window close to E„), see Fig|| Nevertheless, the 
Lyapunov exponent is positive for many initial conditions 
inside the same energy window. 
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FIG. 6. Single particle energy distribution for N — 4, and 
J — 1. Dashed line is obtained dynamically, by integrating 
the equations of motion for a single trajectory with energy 
E* — —1.9745 up to the time t = 10 7 . Full line is the phase 
average distribution obtained for 10 6 different initial condi- 
tions in the energy range [—1.98, —1.97]. 

This lack of ergodicity is also reflected in the lack of 
equipartition. Indeed, defining the average unperturbed 
single particle energy, as 



(5?) = lim 



dt S?(t) 



(59) 



we get that, typically, for trajectories with energy E close 
to E* , neighbors spins do not share the same energy (lack 
of equipartition). In Table I , we show the average ki- 
netic energy per spin for a dynamical trajectory with 
the energy E* (second column), to be compared with 
the average kinetic energy per spin as obtained from the 
phase average distribution (right column) within the en- 
ergy range AE = [-1.98, -1.97]. 

TABLE I. Average kinetic energy per spin. 



spin label i 


(Sf(E*)) 


(S?(AE)) 


1 


-0.7248 


-0.3783 


2 


-0.1387 


-0.3788 


3 


-0.7248 


-0.3788 


4 


-0.1388 


-0.3801 



C. Ergodicity versus dynamical correlations 

In this section we deal with systems having a small 
number of spins (let us say, order of ten). In contrast 
to the case considered in the previous Section, we se- 
lect only those energy values which correspond to both a 
positive Lyapunov exponent and to equipartition among 
different spins; this means, in particular, that "dynami- 
cal" and "statistical" distributions (the former obtained 
by integration of equations of motion for one trajectory, 
the latter by choosing many different initial conditions 
onto the energy surface and performing a phase average) 
are close one to each other. 

Let us focus on two different samples of N = 5 and 
A = 10 interacting spins with strong and chaotic inter- 
action (J = 1). In Figs.0|| we show the SPE-distribution 
for different energies, together with the correspondent 
Gaussian approximations Eq.(f43|). 



o E= -2.294 

— E=-2.294 (Gaussian) 
o E=-4.67 

— E=-4.67 (Gaussian) 
+ E=-3.022 




FIG. 7. Single particle energy distribution n(h\E) for 
N — 5 with interaction J — 1 as a function of single particle 
energy h. Different symbols indicate different energy values 
E, as shown in the window. For few sets of symbols the Gaus- 
sian approximation, as given by Eq.(^3|) has also been shown. 

The first important point is a remarkable deviation 
from the Gaussian approximation, at least for energy 
values E close to the band edge. This simply means 
that the approximations involved in order to get Eq. (E^) 
are no longer valid. One of these approximations was to 
consider a Gaussian shape for the SE in Eq. (0). Ac- 
cording to additional data, the effective SE (statistical 
or dynamical, they are very close one to each other) for 
the same energy value show remarkable deviations from 
the Gaussian. Nonetheless, the substitution of the "true" 
SE in Eq.(fi"^) does not affect considerably n(h\E). This 
amount to say that Eq.([42|) itself does not constitute a 
good approximation for energy values close to the edges. 

Indeed, due to a small number of spins, and to the 
relatively large energy shared by each spin particle, they 
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are strongly correlated and most of the uncorrelation as- 
sumptions (between Hq and V, for example) made in the 
previous chapters are no longer valid. 
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FIG. 8. The same for JV = 10. 

To be more precise, it can be shown that Eq.([l2|) can 
be considered as the "diagonal" approximation of the fol- 
lowing exact relation 



n(h\E) = / dE W E (E ) n E (h\E ) 



where 



n E (h\E ) = P(St = h\H = E ,H = E) 



(60) 



(61) 



In the limit when H and Hq are independent one recov- 
ers Eq.(||) since n E {h\E ) = n (h\E ). The study of the 
correlation kernel Eq.(|6l|) beyond the diagonal approxi- 
mation will be deserved for future investigations. 

On the other side, in the middle of the energy band 
\E\ ~ there is a rough agreement with the Gaussian 
approximation (see Fig.0). Indeed, for such energy val- 
ues, there is no preferred direction of the spin (the energy 
shared by each single spin is relatively small) and Eq. ( Eg ) 
still represents a good approximation. However, the bad 
point is that close to the center of the band the Gaussian 
approximation for n(h\E) is far from an exponential (in 
fact, it is close to a Gaussian!, see Eq.fl3|). Indeed, let 
us remember that one of the conditions in order to ob- 
tain the Boltzman distribution was \E\ ^> 1, which is of 
course not satisfied at the center of the energy spectrum. 

Even if not completely satisfactory from the theoreti- 
cal point of view, but in close analogy with the case of 
N = 100 spins, one could define a temperature as (minus) 
the slope of the fitting straight line to log n(h\E). Oper- 
atively we observe that in the region to the left of their 
intersection point (h ~ —0.6 for N = 5 and h ~ —0.7 for 
N = 10), distributions with different energy value have a 



behavior closer to an exponential. Therefore, it is natu- 
ral to fit the numerical n(h\E) with an exponential only 
to the left of the intersection point. Let us call j3fu the 
inverse temperature obtained in this way. 

Needless to say, when compared with both the statisti- 
cal and microcanonical temperature, one can see impor- 
tant deviations. In Figs.||,|l(] we plot the obtained (3 fa 
as a function of the energy E, together with other defi- 
nitions of temperature for N = 5 and N = 10. 



1.5 



0.5 




FIG. 9. Comparison between different definitions of tem- 
perature for N — 5. Symbols represent the numerically ex- 
tracted (3 fa, while different curves refer to different definitions 
of temperature, microcanonical (full line) and the Gaussian 
approximation (dotted line). Since the latter turns out to be 
dependent on h too, we computed it at ft = —1. 
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FIG. 10. The same for N = 10 spins. 

As one can see, none of the previous definitions seems 
to fit the numerical values. This is not surprising, in 
spite of the fact that only an approximate exponential 
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behavior has been found. We do not have any approxi- 
mate theory able to describe such temperature differences 
when the number of spins is small. Strictly speaking, one 
recognizes the importance of the classical SE in the de- 
scription of the behavior of single particle distribution, 
but it becomes technically complicated to go beyond the 
diagonal approximation, which is correct only when the 
number of particles is sufficiently large and the energy is 
not too close to the "many-body ground state" (bottom 
of the energy spectrum) . 




FIG. 11. f3 fa obtained by fitting the single particle energy 
distribution in the negative part of the single particle spec- 
trum, versus the proposed phenomenological scaling /3 S . Full 
circles represent the numerical data for N = 4, 5, 10, 20, 100, 
and J = 0.1, 0.5, 1 and different energy. Full line is the scaling 
relation. 

The following approximate phenomenological scaling 
relation has been numerically found: 



m 

N 



(1 + a) J 2 
N 



(62) 



where a = 0.7 is the fitting parameter and (E) = 
J dh h n(h\E) . In Fig.[n] for different N, J and energy 
E we show the numerical data and the scaling relation 
Eq.©. 

We have not, for the time being, a theoretical expla- 
nation of this scaling relation. It should be stressed that, 
in the limit TV — > oo and sufficiently large energy (E) the 
second term in the l.h.s. of Eq.(|62|) is negligible with re- 
spect to the first one and Eq.(|5l|) is recovered. We remind 
that it is not possible to take J ^ 1 in this model since 
it becomes integrable, and most of the previous results 
are necessarily wrong. 

The disagreement between the microcanonical temper- 
atures and the statistical one for a small number of parti- 
cles has few important theoretical implications. First of 
all we note that they are obtained in two completely dif- 
ferent ways. While the single particle distribution reflects 
a property of the constant energy surface H = E, the 
microcanonical definition requires a derivative across the 



energy surface. In principle, the knowledge of the correla- 
tion kernel Eq.(|6l|) would solve completely the problem. 
But this in turn requires the knowledge of the (infinite) 
intersections among the H — E surface and the Hq = Ea 
ones. 



VI. CONCLUDING REMARKS 

In this paper we have studied the emergence of Boltz- 
mann's law for the single particle energy distribution 
n(h\E) in a isolated dynamical model of finite number 
of interacting spins. In the limit of a very large num- 
ber of spins this model allows an analytical treatment. 
We have shown that in this strong limit, Boltzmann's 
distribution, indeed, occurs with an effective statistical 
temperature which coincides with that defined by the 
standard microcanonical temperature. The latter is de- 
fined via the total density of states. Since our analytical 
proof is also valid for a strong spin interaction, it is far 
from trivial. Indeed, it is not a priori clear that the SPE 
distribution, which pertains a non-interacting property, 
follows the Boltzmann distribution with the temperature 
determined via the total density of states. 

The above result has been obtained in an approach 
which is very similar to that recently suggested in the 
study of the so-called two-body random interaction model 
[Q-Q. According to this approach, the distribution of 
occupation numbers for single-particles levels in a quan- 
tum many-body system, is directly related to the aver- 
age shape of chaotic eigenstates in the basis of the un- 
perturbed Hamiltonian (the latter can be considered as 
the mean- field part of a system). This means that, in 
fact, there is no needs to know exact eigenstates of huge 
Hamiltonian matrices which takes into account two-body 
interaction between particles. This is due to the chaotic 
nature of eigenstates, which results from the (assumed) 
randomness of the two-body matrix elements. 

An important point of the above approach is that in 
some cases the average shape of eigenstates can be found 
analytically from off-diagonal matrix elements of the to- 
tal Hamiltonian Hq, see Ref. The same happens to 
another important quantity, the strength function, which 
in solid state physics is also known as the local density of 
states (LDOS). These two quantities are related to each 
other, however, so far this relation is not well-studied. 
Numerical data for different models, disordered || and 
dynamical plpl as well, have shown that for not very 
strong interaction, these two quantities are very close one 
to each other. The knowledge of LDOS is very impor- 
tant in many applications. Indeed it gives the informa- 
tion on how energy, initially concentrated in a specific 
unperturbed state, spreads over all other states due to 
the interaction between particles. The inverse width of 
LDOS is, in fact, the effective time of this spread. 
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Until recently, the above two quantities have been dis- 
cussed only in the context of quantum systems. On the 
other hand, in j2^] it was noted that both LDOS and SE 
have a very clear classical analog. The study of the clas- 
sical LDOS and SE have been started in [^D, and the 
first results have shown that the average shape of quan- 
tum eigenstates (and the same for the LDOS) for chaotic 
eigenstates coincides fairly well with the classical coun- 
terpart. For this reason, when studying the occurrence 
of Boltzmann's distribution in our model of interacting 
spins, we have also paid attention to the classical LDOS 
and SE. 

Our results in what concerns classical SE and LDOS 
have shown a nice correspondence to the main findings 
for the quantum model with random two-body interac- 
tion HHH. I n particular, in our classical model we have 
proved the basic expressions for an increase of both to- 
tal energy (^|) and temperature due to statistical 
effects, which in quantum model have been derived di- 
rectly from the shape of exact eigenstates in the unper- 
turbed energy basis. Moreover, new relations ( |32| ) and 
( [40| ) have been found which may be important also in 
quantum systems. 

Numerical data for our model have confirmed main 
theoretical predictions for a large number of particles, 
N = 100. However, it was also found that in spite of a 
relatively large N, one can detect clear deviations for the 
statistical temperature related to the SPE-distribution. 
Detailed analysis shows that these deviations are mainly 
due to dynamical effects of correlations originating from 
a finite number of particles. As was shown, analytical 
results can be interpreted as a "diagonal" approximation, 
that is neglecting correlations between the total and the 
unperturbed Hamiltonian. One should stress that the 
correction to the thermodynamical expressions are not 
only on the order of 1/N as commonly assumed in the 
literature, but depend on the interaction strength as well. 

Drastic deviations from the thermodynamical ap- 
proach have been found for few interacting spins (N = 10 
and N — 5). Namely, both the thermodynamical tem- 
perature defined from the microcanonical relations, and 
the statistical temperature found analytically for a very 
large N, are very different from the approximate temper- 
ature. The latter has been determined numerically from 
the Boltzmann dependence of the SPE distribution, used 
as a fitting expression for low energies of the system. 
A phenomenological expression has been found from the 
analysis of the data, however, without a proper analyt- 
ical explanation. These problems, as well as the devel- 
opments of a semiquantal approach for which quantum 
distribution of occupation numbers is computed with the 
use of classical analog of the shape of eigenstates, will be 
a subject of future investigations. 
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